function [results] = mimo_physical_models(varargin)
% MIMO物理信道模型分析
% 功能：实现MIMO信道物理建模、角度功率谱分析、空间相关性分析
% 输入参数：
%   - carrier_freq: 载波频率 (Hz), 默认 2.4e9
%   - num_paths: 多径数量, 默认 20
%   - angle_spread_tx: 发送端角度扩展 (度), 默认 30
%   - angle_spread_rx: 接收端角度扩展 (度), 默认 60
%   - mean_angle_tx: 平均离开角 (度), 默认 0
%   - mean_angle_rx: 平均到达角 (度), 默认 30
%   - antenna_spacing: 天线间距 (波长), 默认 0.5
%   - Nt: 发送天线数量, 默认 4
%   - Nr: 接收天线数量, 默认 4
% 输出：
%   - results: 结构体，包含所有分析结果和可视化数据
%
% 作者：周勇
% 日期：2024年

%% 参数解析
p = inputParser;
addParameter(p, 'carrier_freq', 2.4e9, @isnumeric);
addParameter(p, 'num_paths', 20, @isnumeric);
addParameter(p, 'angle_spread_tx', 30, @isnumeric);
addParameter(p, 'angle_spread_rx', 60, @isnumeric);
addParameter(p, 'mean_angle_tx', 0, @isnumeric);
addParameter(p, 'mean_angle_rx', 30, @isnumeric);
addParameter(p, 'antenna_spacing', 0.5, @isnumeric);
addParameter(p, 'Nt', 4, @isnumeric);
addParameter(p, 'Nr', 4, @isnumeric);
parse(p, varargin{:});

params = p.Results;

%% 添加路径
addpath('../Common');

%% 获取颜色定义
colors = color_definitions();

fprintf('=== MIMO物理信道模型分析 ===\n');
fprintf('载波频率: %.1f GHz\n', params.carrier_freq/1e9);
fprintf('多径数量: %d\n', params.num_paths);
fprintf('发送端角度扩展: %d 度\n', params.angle_spread_tx);
fprintf('接收端角度扩展: %d 度\n', params.angle_spread_rx);

%% 生成物理信道模型
H_physical = generate_physical_mimo_channel(params.Nt, params.Nr, params.num_paths, ...
    params.angle_spread_tx, params.angle_spread_rx, params.mean_angle_tx, ...
    params.mean_angle_rx, params.antenna_spacing);

%% 分析信道特性
fprintf('\n物理信道模型特性:\n');
fprintf('信道矩阵维度: %dx%d\n', size(H_physical, 1), size(H_physical, 2));
fprintf('Frobenius范数: %.2f\n', norm(H_physical, 'fro'));
fprintf('条件数: %.2f\n', cond(H_physical));

%% 奇异值分解
[U, S, V] = svd(H_physical);
singular_values = diag(S);

%% 角度功率谱计算
angles_tx = -90:5:90;
angles_rx = -90:5:90;
power_spectrum_tx = zeros(size(angles_tx));
power_spectrum_rx = zeros(size(angles_rx));

for i = 1:length(angles_tx)
    power_spectrum_tx(i) = exp(-(angles_tx(i) - params.mean_angle_tx)^2 / (2 * params.angle_spread_tx^2));
    power_spectrum_rx(i) = exp(-(angles_rx(i) - params.mean_angle_rx)^2 / (2 * params.angle_spread_rx^2));
end

%% 空间相关矩阵
R_tx = V * V';
R_rx = U * U';

%% 可视化结果
figure('Name', 'MIMO物理信道模型', 'Position', [100, 100, 1200, 800]);

% 信道幅度响应
subplot(2,3,1);
imagesc(abs(H_physical));
colorbar;
title('信道幅度响应');
xlabel('发送天线');
ylabel('接收天线');

% 信道相位响应
subplot(2,3,2);
imagesc(angle(H_physical));
colorbar;
title('信道相位响应');
xlabel('发送天线');
ylabel('接收天线');

% 奇异值分布
subplot(2,3,3);
plot(singular_values, 'o-', 'LineWidth', 2);
grid on;
title('信道奇异值');
xlabel('模式索引');
ylabel('奇异值');

% 角度功率谱
subplot(2,3,4);
plot(angles_tx, power_spectrum_tx, 'b-', 'LineWidth', 2);
hold on;
plot(angles_rx, power_spectrum_rx, 'r--', 'LineWidth', 2);
grid on;
title('角度功率谱');
xlabel('角度 (度)');
ylabel('相对功率');
legend('发送端', '接收端');

% 发送端相关矩阵
subplot(2,3,5);
imagesc(abs(R_tx));
colorbar;
title('发送端相关矩阵');
xlabel('发送天线');
ylabel('发送天线');

% 接收端相关矩阵
subplot(2,3,6);
imagesc(abs(R_rx));
colorbar;
title('接收端相关矩阵');
xlabel('接收天线');
ylabel('接收天线');

%% 结果打包
results = struct();
results.H_physical = H_physical;
results.singular_values = singular_values;
results.power_spectrum_tx = power_spectrum_tx;
results.power_spectrum_rx = power_spectrum_rx;
results.R_tx = R_tx;
results.R_rx = R_rx;
results.angles_tx = angles_tx;
results.angles_rx = angles_rx;
results.params = params;

fprintf('\nMIMO物理信道模型分析完成！\n');

end

%% 辅助函数
function H = generate_physical_mimo_channel(Nt, Nr, num_paths, angle_spread_tx, angle_spread_rx, mean_angle_tx, mean_angle_rx, spacing)
% 生成基于物理的MIMO信道
H = zeros(Nr, Nt);

% 角度转换为弧度
mean_angle_tx_rad = deg2rad(mean_angle_tx);
mean_angle_rx_rad = deg2rad(mean_angle_rx);
angle_spread_tx_rad = deg2rad(angle_spread_tx);
angle_spread_rx_rad = deg2rad(angle_spread_rx);

for path = 1:num_paths
    % 随机角度
    angle_tx = mean_angle_tx_rad + angle_spread_tx_rad * randn();
    angle_rx = mean_angle_rx_rad + angle_spread_rx_rad * randn();
    
    % 随机复增益
    gain = sqrt(0.5) * (randn() + 1i * randn()) / sqrt(num_paths);
    
    % 发送端阵列响应
    a_tx = exp(-1i * 2 * pi * spacing * (0:Nt-1)' * sin(angle_tx));
    
    % 接收端阵列响应
    a_rx = exp(-1i * 2 * pi * spacing * (0:Nr-1)' * sin(angle_rx));
    
    % 累加路径贡献
    H = H + gain * a_rx * a_tx';
end
end